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Abstract. The Numerical INJection Analysis (NINJA) project is a collabora- 
tive effort between members of the numerical relativity and gravitational wave 
data analysis communities. The purpose of NINJA is to study the sensitivity 
of existing gravitational-wave search and parameter-estimation algorithms using 
numerically generated waveforms, and to foster closer collaboration between the 
numerical relativity and data analysis communities. The first NINJA project 
used only a small number of injections of short numerical-relativity waveforms, 
which limited its ability to draw quantitative conclusions. The goal of the NINJA- 
2 project is to overcome these limitations with long post-Newtonian — numerical 
relativity hybrid waveforms, large numbers of injections, and the use of real detec- 
tor data. We report on the submission requirements for the NINJA-2 project and 
the construction of the waveform catalog. Eight numerical relativity groups have 
contributed 63 hybrid waveforms consisting of a numerical portion modelling the 
late inspiral, merger, and ringdown stitched to a post-Newtonian portion mod- 
elling the early inspiral. We summarize the techniques used by each group in 
constructing their submissions. We also report on the procedures used to validate 
these submissions, including examination in the time and frequency domains and 
comparisons of waveforms from different groups against each other. These proce- 
dures have so far considered only the {(., m) = (2, 2) mode. Based on these studies 
we judge that the hybrid waveforms axe suitable for NINJA-2 studies. We note 
some of the plans for these investigations. 



1. Introduction 

A new generation of laser interferometric gravitational-wave detectors (Advanced 
LIGO [1-3], Advanced Virgo [4,5], and LCGT [6]) is presently under construction. 
These second-generation detectors will have an order of magnitude increase in 
sensitivity over first generation instruments and will be sensitive to a broader range 
of gravitational-wave frequencies. One of the most widely anticipated sources for 
this global network of observatories is the inspiral, merger and ringdown of a binary 
containing two black holes [7]. Detection of such a binary black hole coalescence will 
allow astronomers and astrophysicists to directly observe the physics of black-hole 
spacetimes and to explore the strong-field conditions of Einstein's theory of general 
relativity [8]. 

The ability of gravitational-wave astronomers to use the new generation of 

observatories to detect and study binary black hole coalescence depends on the quality 
of search and source-parameter measurement algorithms. These algorithms rely on 
the physical accuracy of the underlying theoretical waveform models. Developing 
and testing the algorithms required to achieve the goals of gravitational-wave 
astronomy demands close interaction between the source-modeling and data-analysis 
communities. The Numerical INJection Analysis (NINJA) project was created in 
2008 to bring these communities together and to use the recent advances in numerical 
relativity (NR) [9] to test analysis pipelines by adding physically realistic signals 
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to detector noise in software. We describe such additions of signals into noise as 
"injections." 

The first NINJA project (NINJA-1) [10] considered a total of 23 numerical 
waveforms, which were injected into Gaussian noise colored with the frequency 
sensitivity of first-generation detectors. These data were analyzed by nine data- 
analysis groups using both search and parameter-estimation algorithms. However, 
there were two major limitations to the NINJA-1 analysis: First, to encourage 
broad participation, no length or accuracy requirements were placed on the numerical 
waveforms. Consequently, many of these waveforms were too short to inject over an 
astrophysically interesting mass range without introducing artifacts into the data. The 
lowest mass binary considered in NINJA-1 had a total mass of 35M0, whereas the 
mass of black holes could extend below SM© [11, 12]. In NINJA-1, the waveforms were 
only inspected for obvious, pathological errors and no cross-checks were performed 
between the submitted waveforms; and, therefore it was difficult to assess the physical 
fidelity of the results. Second, the NINJA-1 data set contained stationary noise with 
the simulated signals already injected into the data. The data set contained only 
126 simulated signals, which precluded detailed statistical studies of the effectiveness 
of search and parameter estimation algorithms. Finally, since the data set lacked 
the non-Gaussian noise transients present in real detector data, it was not possible 
to fully explore the response of the algorithms in a real search scenario. Despite 
these limitations, NINJA-1 successfully removed a number of barriers to collaboration 
between the source-modelling and data-analysis communities and demonstrated where 
further work is needed. The goal of the second NINJA project (NINJA-2) is to address 
the deficiencies of NINJA-1 and to perform a systematic test of the efficacy of data- 
analysis pipelines in real-world situations in preparation for Advanced LIGO and 
Virgo. 

This paper reports on the improvements we have made to the NINJA analysis to 
address the first of the limitations described above — the accuracy of the numerical 
waveforms. We present the NINJA-2 waveform catalog and describe the results of the 
procedures we have used to validate these data. NINJA-2 places requirements on the 
accuracy and length of the contributed waveforms and we have performed systematic 
cross-checks of the submitted waveforms. Each binary black hole simulation in the 
NINJA-2 catalog must include at least five orbits of usable data before merger, i.e., 
neglecting the initial burst of junk radiation. The NR waveform amplitude should be 
accurate to within 5%, and the phase (as a fimction of gravitational- wave frequency) 
should have an accumulated uncertainty over the entire inspiral, merger and ringdown 
(of the numerical simulation), of no more than 0.5 rad. We also require that numerical 
simulations arc "hybridized" to post-Newtonian (pN) waveforms so that the resulting 
waveforms contain enough cycles to allow injections at M > lOM© in early Advanced 
LIGO data. The continued advances in numerical simulations have also allowed us 
to study a somewhat larger region of the signal parameter space; however we have 
restricted our attention to non-precessing binaries for reasons discussed in Sec. 2 below. 

A subsequent paper will describe the results of using the NINJA-2 waveforms 
to study search and parameter estimation algorithms in real detector data. Since 
data from the second-generation detectors is not yet available, the NINJA-2 analysis 
will use data from the first-generation detectors re-colored so that it has the frequency 
response expected in the first observing runs of the Advanced LIGO (aLIGO) detectors 
— referred to as "early aLIGO" [13]. Similar noise curves will be used to simulate 
the Advanced Virgo detector. NINJA-2 analyses will use these noise models to ensure 
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Figure 1. The hybrid q = m\/m2 = 2, non-spinning MayaKranc waveform 
scaled to various total masses shown against the early and zcro-detuncd-high- 
power aLIGO noise curves, shown as amplitude spectral densities, the square root 
of the power spectral densities. The triangles represent the starting and ending 
frequencies of the post-Newtonian hybridization region, given in table 1. The 
total mass of the binary is scaled so that the hybridization region ends at 100 Hz, 
40 Hz, and 10 Hz. The amplitude of the signal is scaled so that it represents an 
optimally oriented binary at a distance of 1 Gpc from the detector. The early 
aLIGO sensitivity is used to compute the signal-to-noise ratio p. 



that existing algorithms are optimal when second-generation detectors come online in 
2015. The results in this paper use the early aLIGO sensitivity curve (cf. Fig. 1) to 
study the accuracy of the submitted waveforms. The ultimate sensitivity of Advanced 
LIGO is expected to be significantly better than this curve. To allow the waveforms to 
be used in studies using more sensitive noise curves, we have also performed accuracy 
studies using the aLIGO zero- detuned, high-power [14] sensitivity curve. Fig. 1 
shows the two aLIGO sensitivity curves, characterized by their Amplitude Spectral 
Densities (ASD) overlaid with one of the contributed NINJA-2 waveforms. This 
figure demonstrates that hybridization is necessary to allow scaling of the numerical 
waveforms to astrophysically interesting masses, and a portion of the present paper 
studies the hybridization methods used to construct the NINJA-2 waveforms. 

The remainder of this paper is organized as follows: Section 2 describes in more 
detail the accuracy requirements that we have placed on NINJA-2 simulations and 
presents an overview of the waveform catalog showing the regions of the binary black 
hole parameter space covered. Section 3 gives an overview of the numerical methods 
used to construct the numerical relativity waveforms and Sec. 4 describes the methods 
that we have used to hybridize the numerical simulations to pN waveforms. The 
pN waveforms themselves are summarized in the Appendix. Section 5 describes the 
methods and results of the comparisons we have performed between the waveforms. 
Based on these comparisons, we judge the hybrid waveforms suitable for the NINJA- 
2 project. Section 6 summarizes our findings and suggests directions for future 



NINJA-2 waveform catalog 



5 



improvements of the catalog, in particular the study of higher order modes in the 
waveforms. 

2. Overview of the Waveform Catalog 

Binary black holes formed from the evolution of massive stars are expected to have 
circularized before their gravitational-wave frequency reaches the sensitive band of 
ground-based detectors, and so we only consider circular (non-eccentric) binaries. In 
the NINJA-2 project we further restrict our attention to binaries that do not undergo 
precession, i.e., where the spins of the black holes either vanish or are parallel (or 
anti-parallel) to the binary's orbital angular momentum. We do this for two reasons: 
(i) in trying to understand the complex phenomenology of the binary parameter 
space, we prefer to tackle first a simpler subset, which nonetheless captures the main 
features of binary inspiral and merger; (ii) the processing-binary parameter space has 
been sampled by only a handful of numerical simulations. The numerical-relativity 
community is currently exploring the space of prcccssing binaries, for example through 
the numerical relativity-analytical relativity (NRAR) project [15]. Such waveforms 
will be used in future NINJA projects that explore precession. 

The parameters of the black- hole binaries we consider in NINJA-2 are the mass of 
each black hole, mi and m2, or equivalently the total mass M = toi+TO2 and mass ratio 
g = TOi/m2, and the dimensionless spin- magnitude of each black hole, Xi = Sx/m\ 
and X2 = S2/m2- The total mass sets the overall scale of the system, and can be 
factored out to leave a three-dimensional parameter space, {q, xi, X2\- Figure 2 shows 
the coverage of parameter space (details in Sec. 3). The sampling is coarse; while the 
waveforms will provide invaluable information within the NINJA-2 project, we expect 
that ultimately a more uniform coverage of parameter-space by a much larger number 
of configurations will be necessary. 

The NINJA-2 requirement of five pre-merger orbits is at the low end of 
estimates of sufficient waveform lengths for the construction of accurate hybrid PN- 
NR waveforms, as discussed in [16-20], but we expect these to be acceptable for 
the goals of the NINJA-2 project. The 5% amplitude and 0.5 rad phase accuracy 
requirements were formulated with typical current waveforms in mind, for example 
those studied in the Samurai project [21] and studies performed in preparation for 
the NR-AR collaboration project. These requirements are consistent for waveforms 
of similar lengths but may not be directly applicable to much longer waveforms. For 
example, in the 25-orbit SpEC simulations with dimensionless spins Xi = 0.97 the 
highest- and second-highest-resolution data differ by roughly 0.6 rad at merger [22]. 
Note that because this phase-error accumulates over 20 additional inspiral orbits, this 
waveform would easily satisfy the NINJA-2 phase requirement if it were truncated to 
minimally meet the NINJA-2 length requirement (although such a truncation would 
decrease the accuracy of the hybridized waveform). 

We require that the hybridized waveforms in the NINJA-2 catalog are long enough 
to span the sensitivity bands of the advanced LIGO and Virgo detectors in their 
early operation. Specifically, when rescaled to 10 the hybrids must begin with 
a gravitational wave frequency of 20 Hz or lower, i.e. a starting GW frequency of 
Mbj < 0.006. This requires extending the NR waveforms to lower frequencies (i.e. 
more inspiral cycles) by attaching a pN inspiral waveform onto the early portion 
of the NR waveform to produce a hybrid pN-NR waveform. We require that the 
hybridization be performed at a gravitational-wave frequency of Mw22 < 0.075, where 
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Figure 2. Mass ratio q and dimensionless spins Xi of the NINJA-2 hybrid 
waveform submissions. 

Muj22 is the frequency of the {i,m) — (2, ±2) harmonic. In practice, hybridization 
fits were performed over a frequency range as summarized in Sec. 4, and the average 
frequency, and frequency of the average time of the fitting interval were always chosen 
below Mu}22 < 0.075, with two exceptions as seen in Table 1: The nonspinning Llama 
waveforms at mass ratios g = 1, 2. As seen in Fig. (7) these do however show excellent 
overlaps with comparison waveforms. 

The waveforms were submitted with the complex GW strain function — ih^ 
decomposed into modes using spin-weighted spherical harmonics "'^Yim of weight 
s = —2. Although most of the power is in the {£, m) = (2, ±2) modes, we encouraged 
(but did not require), the submission of additional subdominant modes. The accuracy 
studies in this paper focus on the (^, m) = (2, 2) mode; further work is required to 
study the accuracy of the contributed subdominant modes. A total of 63 waveforms 
from 8 groups were contributed to the NINJA-2 catalog. There are 46 distinct 
numerical waveforms; some of these waveforms have been hybridized with multiple 
pN waveforms. The NINJA-2 catalog is summarized in Table 1, and a map of the 
parameter values is shown in Fig. 2. In the next section, we describe in more detail 
the numerical methods used to generate these waveforms and present additional plots 
in Figs. 3 and 4. 

3. Numerical Methods 

3.1. Summary of contributions 

The NINJA-2 data set contains both hybrid and original numerical relativity 
waveforms, in a data format that is summarized in Sec. 3.2 below, and described 
in detail in Ref. [60]. The contributed waveforms cover 29 different black hole 
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q 


Xi 


X2 


Submission 


lOOOe 


100 Mw 


#NR 


pN 












hyb. range 


cycles 


Approx 


1.0 


-0.95 


-0.95 


SpEC [22- 31] 


1.00 


3.3 - 


4.1 


18.42 


Tl 


1.0 


-0.85 


-0.85 


BAM [32 37] 


2.50 


4.1 - 


4.7 


12.09 


T1,T4 


1.0 


-0.75 


-0.75 


BAM [32-37] 


1.60 


4.1 


4.7 


13.42 


T1,T4 


1.0 


-0.50 


-0.50 


BAM [32-37] 


2.90 


4.3 - 


4.7 


15.12 


T1,T4 


1.0 


-0.44 


-0.44 


SpEC [23,24,27-30,38] 


0.04 


4.3- 


5.3 


13.47 


T4 


1.0 


-0.40 


-0.40 


Llama [39-41] 




6.1 - 


8.0 


6.42 


T1,T4 


1.0 


-0.25 


-0.25 


BAM [32-37] 


2.50 


4.5 - 


5.0 


15.15 


T1,T4 


1.0 


-0.20 


-0.20 


Llama [39-41] 




5.7- 


7.8 


8.16 


T1,T4 


1.0 


0.00 


0.00 


BAM [32-37] 


1.80 


4.6 - 


5.1 


15.72 


T1,T4 








GATech [42-48] 


3.00 


5.5 - 


7.5 


9.77 


T4 








Llama [40, 41] 




5.7- 


9.4 


8.30 


F2 








SpEC [23, 24, 27-30] 


0.05 


3.6 


4.5 


22.98 


T4 


1.0 


0.20 


0.20 


GATech [42-48] 


10.00 


6.0 - 


7.5 


10.96 


T4 


1.0 


0.25 


0.25 


BAM [32-37] 


6.10 


4.6- 


5.0 


18.00 


T1,T4 


1.0 


0.40 


0.40 


GATech [42-48] 


10.00 


5.9 - 


7.5 


12.31 


T4 








Llama [39-41] 




7.8 - 


8.6 


6.54 


T1,T4 


1.0 


0.44 


0.44 


SpEC [23, 24, 27-30, 38] 


0.02 


4.1 - 


5.0 


22.39 


T4 


1.0 


0.50 


0.50 


BAM [32-37] 


6.10 


5.2 - 


5.9 


15.71 


T1,T4 


1.0 


0.60 


0.60 


GATech [42-48] 


12.00 


6.0 - 


7.5 


13.63 


T4 


1.0 


0.75 


0.75 


BAM [32-37] 


6.00 


6.0- 


7.0 


14.03 


T1,T4 


1.0 


0.80 


0.00 


GATech [42-48] 


13.00 


5.5 


7.5 


12.26 


T4 


1.0 


0.80 


0.80 


GATech [42-48] 


6.70 


5.5 - 


7.5 


15.05 


T4 


1.0 


0.85 


0.85 


BAM [32-37] 


5.00 


5.9- 


6.9 


15.36 


T1,T4 








UIUC [49] 


20.00 


5.9 


7.0 


15.02 


Tl 


1.0 


0.90 


0.90 


GATech [42-48] 


3.00 


5.8 - 


7.5 


15.05 


T4 


1.0 


0.97 


0.97 


SpEC [22-25,27-31] 


0.60 


3.2 - 


4.3 


38.40 


T4 


2.0 


0.00 


0.00 


BAM [32-37] 


2.30 


6.3 - 


7.8 


8.31 


T1,T4 








GATech [42-48] 


2.50 


5.5 - 


7.5 


10.42 


T4 








Llama [40, 41] 




6.3- 


9.4 


7.47 


F2 








SpEC [23-25,27-31,50] 


0.03 


3.8 


4.7 


22.34 


T2 


2.0 


0.20 


0.20 


GATech [42-48] 


10.00 


5.6 - 


7.5 


11.50 


T4 


2.0 


0.25 


0.00 


BAM [32, 36] 


2.00 


5.0- 


5.6 


15.93 


T1,T4 


3.0 


0.00 


0.00 


BAM [32-37] 


1.60 


6.0 - 


7.1 


10.61 


T1,T4 








SpEC [23-25,27-31,50] 


0.02 


4.1 - 


5.2 


21.80 


T2 


3.0 


0.60 


0.40 


FAU [51-54] 


1.00 


5.0- 


5.6 


18.89 


T4 


4.0 


0.00 


0.00 


BAM [32-37] 


2.60 


5.9 - 


6.8 


12.38 


T1,T4 








LEAN [55, 56] 


5.00 


5.1 - 


5.5 


17.33 


Tl 








SpEC [23-25,27-31,50] 


0.03 


4.4- 


5.5 


21.67 


T2 


6.0 


0.00 


0.00 


SpEC [23-25,27-31,50] 


0.04 


4.1 


4.6 


33.77 


Tl 


10.0 


0.00 


0.00 


RIT [57-59] 


0.40 


7.3 - 


7.4 


14.44 


T4 



Table 1. Summary of the NINJA-2 waveform catalog. Given are mass-ratio 
q = mi/m2, magnitude of the dimcnsionless spins Xi = Si/rnf, numerical code, 
orbital eccentricity e, frequency range of hybridization in Muj, the number of 
numerical cycles from the middle of the hybridization region through the peak 
amplitude, and the post-Newtonian Taylor-approximant(s) used for hybridization. 
All pN approximants include terms up to 3.5 pN-order, see the Appendix. 
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configurations modeling low-eccentricity inspiral, the mass ratio q = m\/m2 ranges 
from 1 to 10, and the simulations cover a range of non-prcccssing spin configurations. 

The initial frequency uj of the m) = (2,2) mode for the numerical waveforms 
ranges from 0.035/M to 0.078/M, where M denotes the sum of the initial black- 
hole masses, with both mean and median values of 0.048/Af. Table 1 lists a few key 
parameters that distinguish the waveforms, and introduces short tags for the different 
contributors: 

(i) Two groups use versions of the BAM code, "BAM" labels the Cardiff-Jena-Palma- 
Vienna collaboration [32,33,36,37,54], and "FAU" the contribution from the 
Florida Atlantic group [32, 36, 54, 61]. 

(ii) LazEv is the RIT code [58, 62, 63]. 

(iii) LEAN has been developed by Ulrich Sperhake [55, 56] . 

(iv) Two contributions use the Llama code [40,41]. Llama- AEI is the contribution of 
the AEI group [40,41], LLama-PC is the Palma-Caltech contribution [39]. 

(v) GATech is the Georgia Tech group, using the MayaKranc code [48,64]. 

(vi) SpEC for the Cornell-Caltech-CITA collaboration code [24, 27-29], 

(vii) UIUC stands for the University of Illinois at Urbana-Champaign team [65-67]. 

The numerical codes follow either of two approaches to solving the Einstein 
equations (see [68] for a review). The SpEC code employs the generalized harmonic 
formulation (see e.g. Refs. [69,70]) with gauge conditions adapted to for black hole 
binaries [25,27,29]. SpEC employs black hole excision to remove singularities in 
the interiors of the black holes from the computational domain. SpEC's initial-data 
(also using black hole excision) [71-73] is constructed with a pseudo-spectral elliptic 
solver [23]. 

All other codes use the BSSNOK formulation of the Einstein evolution 
equations [74-76] with hyperbolic evolution equations for the lapse and shift in 
the moving punctures formalism [58,77]. The 1+log slicing condition for the lapse 
function [78] is "singularity-avoiding": the time slices freeze in before reaching the 
singularity in the black hole. This makes it possible to avoid the use of black hole 
excision techniques [79 82], when evolving the shift vector field /3* according to the 
F-driver condition [83, 84] (extended to the moving puncture approach which allows 
for some free parameters which groups tune individually). 

A significant amount of computational infrastructure is shared between a number 
of codes. With the exception of the two groups using BAM, all other moving-puncture 
codes are based on the Cactus computational toolkit [85,86], the Carpet mesh- 
refinement code [87,88] or the EinsteinToolkit infrastructure [89,90]. The Cactus- 
based codes also use the same apparent horizon finder code (AHFinderDirect) [91]. 
The codes Llama, LazEv, Lean and MayaKranc all use the same pseudospectral solver 
for the Einstein constraint equations [92], and BAM \iscs a variant thereof [32]. 

We will only very briefly summarize the main features of the numerical methods 
and codes, as such information is generally available elsewhere (see references above 
and the NINJA-1 paper [10]). There are two important exceptions: Two groups use 
the new Llama code, which is based on a multipatch decomposition of the numerical 
grid, and uses spherical coordinates in the outer zones of the grid, similar to the 
SpEC code. This allows a more efficient treatment of the wave zone, and to causally 
disconnect the outer boundaries from the wave extraction; in addition Llama uses 
characteristic extraction to extract the waves directly at null infinity [93,94]. The 
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other important new development is the inclusion of two configurations with BH Kerr 
parameters 0.95 and 0.97 [22, 26], which was made possible by using superposed Kerr- 
Schild initial data [73,95] in contributions from the SpEC collaboration. Suitably 
conformally curved initial data allows BH spins beyond the Kerr parameter of w 0.93, 
which is the maximum attainable using conformally flat initial data ([73] and the 
references therein) . For further details on all the numerical codes used, we refer to the 
code references listed above, and the recent overview papers on the numerical solution 
of the binary black hole problem [68, 96, 97] . 

3.2. Data format of contributions 

All contributions followed the format specified in Ref. [60]. The data format consists 
of metadata files and data files with spherical harmonic modes. The metadata 
specify the physical parameters of the BH binaries, such as mass ratio and spins, the 
initial frequency of the {d., m) — (2, 2) spherical-harmonic mode, eccentricity, and also 
authors, bibliographical references, as well as mimerial methods used. This metadata 
format has been significantly extended since the first NINJA project to contain more 
information about the numerical simulations. For NINJA-1, the waveform data were 
stored as 3-column ASCH tables, listing the time at equidistant steps, and real and 
imaginary parts of the strain. For the long hybrid waveforms in NINJA-2, this format 
is not efficient; rather we store the time, amplitude and phase of the modes. The 
amplitude and phase as functions of time exhibit much less temporal structure than the 
complex waveform's oscillatory behavior. Therefore, amplitude and phase at arbitrary 
time can easily be recovered by interpolation from a drastically reduced number of 
time steps. Consequently, data may be provided with unequal time-spacing with only 
as many steps as required to accurately regenerate the original hybrid waveform with 
simple linear interpolation. 

3.3. Initial data and eccentricity 

Specifying initial data for black hole binaries in a non-eccentric inspiral is by itself 
a non-trivial problem. The elliptic constraint equations of general relativity need to 
be solved numerically, and the free data, which serve as input to these equations and 
select a specific configuration of black holes, have to be chosen in a judicious way (for 
a general overview see, e.g., Ref. [98]). 

The moving puncture and generalized harmonic codes differ in the way they 
specify the free data for the constraint equations, and correspondingly in how they 
encode the black hole parameters. The codes based on the "moving puncture" 
approach use puncture initial data [99-102] to model black holes, resulting in initial 
data that contain a separate asymptotically flat end within each black hole. The 
lapse and shift fields, which determine the coordinate gauge, are initially set to 
trivial values, and quickly pick up values that keep the geometry of the black holes 
smooth and almost time independent. The SpEC code uses quasi-equilibrium excision 
initial data where the interior of the black-hole horizons has been excised from the 
numerical grid [71, 72, 103]. The constraint equations are solved using the conformal- 
thin-sandwich formulation of the initial- value problem [104,105]. These data also 
specify an initial lapse and shift, thus the evolutions can already be started in an 
appropriate coordinate gauge. 

All of the codes solve the elliptic constraints with pseudo-spectral numerical 
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codes [23,92], sharing numerical infrastructure as noted above. Initial data always 
correspond to the ccntcr-of-mass rest frame, such that the net linear momentum 
vanishes initially. All codes take input parameters that determine the individual 
black- hole masses rrij, spins Si, momenta Pi and coordinate separation D of the 
black holes. We note however that in the dynamical strong field regime of general 
relativity, definitions of mass, spin, and linear momentum of individual black holes 
are ambiguous. In our case, initial separations of the black holes are large enough to 
match with post-Newtonian waveforms to construct hybrids, and such ambiguities are 
in some sense hidden in this matching procedure and associated errors. 

In order to achieve non-eccentric inspiral, appropriate input parameters have to 
be found, i.e., appropriate momenta have to be chosen for given masses, spins, and 
separation. No precise definition of eccentricity is available in general relativity, and 
one therefore usually resorts to quantities inspired from Newtonian gravity, which 
quantify those oscillations in the black hole tracks or wave signal that are associated 
with eccentricity, see e.g., the discussion in Ref. [106]. The initial parameters are 
determined by a number of different methods, using either an initial guess from 
a standard post-Newtonian or efFcctive-one-body (EOB) approximation based on 
Refs. [107-109], or adding an iterative procedure to further reduce the eccentricity, 
based on Refs. [29,50, 110, 111]. Measured orbital eccentricities are listed in Table 1, 
and have been checked for consistency by an independent estimate of the eccentricity 
based on the submitted GW information. 

3.4- Gravitational-wave extraction 

The gravitational-wave signal can only be defined unambiguously at null infinity. The 
first code that is capable of computing the wave signal at null infinity for black 
hole coalescences has become available since the NINJA-1 project, and is based 
on combining a characteristic extraction method with the Llama code [40,41], and 
several waveforms have been contributed using this code. Also, these calculations have 
provided rigorous error estimates for procedures where the wave signal is extracted 
at finite radius, or at a sequence of radii, combined with an extrapolation to infinite 
extraction radius [112], thus providing justification to such approximations. 

Information about extraction radius were included for most contributions. The 
extraction radii ranged from 75M to 500 Af, with a median of 90M. For four 
contributions the GW signal was read off at null infinity using characteristic 
extraction [40, 113], roughly ten contributions extrapolated the GW signal to infinity 
from a number of extraction radii. Two techniques are used to extract an approximate 
gravitational waveform at finite extraction radius: SpEC uses the Regge-Wheeler- 
Zerilli method [112, 114-117] to compute the strain directly; all other contributions 
use the Newman-Penrose curvature scalar ^4 [118]. Computation of the strain from 
•^4 requires two time integrations, which requires the proper choice of the constants 
of integration, and may require further "cleaning procedures" to get rid of artifacts 
resulting from the finite extraction radii, as discussed in detail in Ref. [119]. Several 
techniques were used in practice, based on time domain [120], or frequency domain 
methods such as [119], or heuristic methods of suppressing low frequency Fourier 
modes, or combining the latter with the method of [120] as in the BAM submissions. 
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3. 5. Numerical Methods, accuracy, and sources of error 

The numerical rtlgoiitlims employed in the various codes agree in many details: For 
the time discretization, all codes use fourth or fifth order accurate explicit algorithms 
based on the method of lines. The consistent experience of different groups is that 
finite difference errors are dominated by the spatial discretization, correspondingly all 
the moving puncture codes use at least sixth order finite differencing for the spatial 
discretization. The SpEC code uses a multi-domain pseudospectral method with a 
large number of domains, which shows exponential convergence. 

The moving puncture codes use a simple hierarchy of fixed refinement boxes 
which follow the motion of the black holes, and information between the boxes is 
communicated using buffer zones and a variant of Berger-Oliger mesh-refinement. 
Interpolation between refinement levels is performed with spatial polynomial 
interpolation of fifth (Llama, UIUC, LazEv, Lean and MayeiKraiic) or sixth order (BAM). 
Time interpolation at mesh-refinement boundaries introduces second-order errors, 
which does however not appear to dominate numerical errors. 

4. PN-NR Hybrid Waveforms 

While post-Newtonian (pN) methods accurately approximate GW signals throughout 
early inspiral, they become increasingly unreliable towards late inspiral and 
merger [121]. Numerical relativity (NR) simulations are capable of accurately 
computing inspiral, merger, and ringdown portions of binary coalescence [58, 77, 96, 
122], but are too computationally expensive to extend far into the inspiral regime [27]. 
Hybrid waveforms arc the result of smoothly blending together pN and NR waveforms 
to form a waveform that covers the full BBH dynamics. For the NINJA-2 project, 
each NR group has produced their own hybrid waveform after ensuring that the pN 
portions all agree. We expect some systematic errors resulting from errors in the pN 
approximation, the choice of blending region, and the hybridization method [16, 18- 
20,27,123,124]. The NINJA-2 hybrid waveforms do not contain effective-one-body 
extensions of pN approximants, although these have also been used to model complete 
waveforms (e.g. [19,125,126]). 

Hybridization uses least-squares fits to determine the extrinsic parameters for the 
pN waveform [124, 127, 128]. In general, this is accomplished by evaluating 



where T represents waveform data relating to strain [e.g., h{t) = h+{t) — ihx{t), 
aTg[h{t)] or h{f)]. If T is derived from the time domain, then s = if T is in the 
frequency domain, then s = f. For either case, [.si,S2], chosen within the domain of 
both the pN and NR data sets, defines the integration interval and, in most cases, the 
blending region. The vector u denotes the set of pN-parameters over which the fitting 
is performed. For example, u = (tshiit, ^sinft, m) corresponds to adjusting time- and 
phase- shift and the mass ratio of the pN waveform to match the NR waveform. The 
best-fit parameters are denoted by u*. The amplitude scaling factor, a, is often fixed 
to a = 1, but may be included in the fitting parameters [127]. Finally, in the limit 
Si — >■ S2, this procedure reduces to enforcing equality of TpM and Tnr at si = S2, as 
well as equality of the first derivative. 

More explicitly, hybridization may be performed via the following algorithm: 




(1) 
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(i) Choose [si,S2] within the pN and NR data sets. Ideally, [si,S2] is sufficiently 
early so that both pN and NR sets should be accurate. 

(ii) Evaluate Eq. (1); apply {u* ,a*} to the pN data set, resulting in T*f^. Measure 
error quantities relating to fit. 

(iii) If desired, adjust [si, S2] and iterate (i) and (ii), to find a preferred interval [s^, S2]. 

(iv) Defining a monotonic function z{s) such that z{s < = and z{s > Sb) = 1, 
the hybrid is given by 

THyb(s) = [1 - z{s)] + a* z{s) Tnr . (2) 
Note that the transition region [sa,Sb\ is generally taken to be a sub-interval of 
[sj, S2], sometimes consisting of a single point. 

This general formalism leaves many decisions open, and the following specific 
choices were made during hybridization: 

• For the SpEC waveforms the integrand of Eq. (2) was taken to be the square of the 
phase-difference only [128], T = SLTg[h(t)], and maximation was performed over 
tshih and (/ishift only, without any adjustments to the amplitude (a* = 1). The 
time-interval [ti,t2] was chosen to correspond to the frequencies listed in Table 1, 
without any iterative adjustments of ti and ^2- 

• The RIT waveforms similarly use T = a.ig[h{t)] for Eq. (1), but employ the limit 
as ti t2 to determine u = (ishift) ^^shift) at Moo = 0.075. Then, the transition 
function is znnit) = x(tf [6x{tf - 15x{t) + 10], where x{t) := (t - ti)/(i2 - h) 
for a finite interval {ti ^ ^2), guaranteeing behavior at t = ti and ^2 [129]. 

• For the Lean waveforms hybridization is performed in a similar way [130], using 
the transition function z{t) = 70a;'' — 315.T*^ + 540a;^— 420a;^ + 126.T''', and individual 
mode amplitudes of the PN waveforms are rescaled such that their average over 
the matching window agrees with the numerical result. 

• The UIUC waveforms use T = h{t), a = 1, using as free parameters the initial 
PN phase and orbital angular frequency. Equation (2) is then used to construct 
the final hybrid, with z{t) — {t — ti) / (t^ — ti) as in [131]. The time-interval [^1,^2] 
was chosen to correspond to the UIUC frequencies listed in Table 1, without any 
iterative adjustments of ti and ^2- 

• The GATech hybridization follows [128] and is done in the time domain with 
T = h{t) and u = (tshift) ^^shift)- Equation (1) is evaluated over {u, a} and then 
equation (2) is used to construct the final hybrid, with z{t) = {t — ti)/{t2 — ti). 
The fitting intervals are given in Table 1. 

Figure 3 and 4 show exemplary plots of the resulting hybrid waveforms. For 
aligned spins, orbital hangup extends the inspiral to smaller separation and higher 
frequencies. This is apparent in Fig. 3 in that the last ten GW cycles (as indicated by 
the small circles) take less time, and in Fig. 4 by the shift toward higher frequencies. 

5. Validation and compeirison of hybrid waveforms 

Each NR group verified that their waveforms met the minimum NINJA-2 requirements 
before submission, as described in Sec. 1. Once submitted, a series of checks were 
performed in order to validate the waveforms against each other. In the first stage 
the post-Newtonian expressions and codes were compared against each other and the 
literature. This resulted in a set of codes in various languages producing waveforms 
that agree well in both phase and amplitude (see Appendix A). 
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Figure 3. Summary of all submitted hybridized waveforms (r/M)h+ The 

X-axis shows time in units of M and the y-axis shows the real part of the {£, m) = 
(2,2) component of the dimensionless wave strain {r/M)h = {r/M){h^ — ihx)- 
The top group shows equal-mass equal-spin waveforms. The middle group shows 
unequal-mass and zero-spin waveforms, and the bottom group show unequal spin 
waveforms. The black circles indicate 10 and 20 GW cycles measured from the 
waveform peak. The hybridization frequency range is shown in the yellow line. 
The post-Newtonian part is shown in red line and the NR portion occurring after 
hybridization is shown in blue line. 
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Figure 4. Sample frequency domain plot Shown are the plots of \h(f)\\/J 
of (l, m) = (2, 2) mode for all equal-mass equal-spin waveforms. The waveforms 
are scaled to 10 Mq and axe placed at 100 Mpc. The Fourier transforms show 
monotonic behavior in the spin parameter x = {iXl +X2)/(? + l) highlighting the 
orbital hang-up effect due to spin. The vertical line indicates 20 Hz, the required 
upper bound on the initial frequency of the hybrids. 

5.1. Time- domain and frequency- domain checks 

In the second stage of validation we examined the {£, m) = (2, 2) mode of the hybrid 
waveforms. We first plotted the last 40 cycles of each waveform — enough to include 
the full NR portion, the hybridization region, and some of the pN portion - and looked 
for any anomalies such as "kinks" caused by the hybridization procedures. Similar 
visual checks were performed on the amplitudes of the Fourier transforms of the entire 
waveforms. This process identified a few issues, including a bug in one hybridization 
code, which were then corrected. While this was useful, it was only possible due 
to the relatively small number of waveforms in NINJA-2 and the fact that only the 
, m) = (2; 2) mode was examined. For future NINJA projects it will be necessary 
to automate this process, the NINJA-2 data analysis may suggest methods for such 
automation. 



5.2. Overlap Comparisons 

In this check the waveforms were compared against each other using matched-filtering 
techniques. The inner product between two real waveforms s\{t) and S2{t) is defined 
as 

(.,|s,)=4jR^°°d/^i^Hp (3) 
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Figure 5. Left: Overlap time series between tlie GATecli and SpEC equal-mass, 
non-spinning waveforms scaled to 20Mq at different sample rates. At 4096 Hz 
the sample point nearest the maximum is sufficiently far that the overlap is 
significantly underestimated, where we are interested in differences to one part in 
10~^. Right: Overlaps between the same waveforms as a function of mass, at 
different sample rates. All overlaps are calculated against the early aLIGO noise 
curve. 



where S,i{.f) is the power spectral density (see Fig. 1). The overlap is then obtained 
by normalization and maximization over relative time and phase shifts, At and A0. 

(si I S2) := max , ^'"^ ' '"'^ ^ . (4) 

There is an important subtlety involved in calculating these overlaps. Usually, all 

possible time-shifts are evaluated simultaneously by a suitable combination of fast 
Fourier transforms. This results in overlaps at time-shifts at discrete values of At 
spaced by the inverse sampling frequency. The result of time-maximization is then 
taken to be the maximum of this discrete time-series. When comparing two very 
similar waveforms, such as two NR simulations of the same system, the overlap 
function becomes very sharply peaked. In units where G = c=l, IMq = 4.93x 10~®s, 
and in these units time shifts of well below IM can lead to significant changes in the 
overlap. It is therefore imperative that the sample rate used in calculating the overlap 
be large enough to find the true maximum. 

This issue is demonstrated in Fig. 5. The left panel shows the overlap function 
resulting from the comparison of two equal-mass, non-spinning waveforms sampled at 
four diff'erent rates. At 4096 Hz the peak of the function is missed and the overlap is 
underestimated. The consequence of this is illustrated on the right of Fig. 5 which 
shows that the overlap as a function of mass exhibits oscillations. As the mass changes, 
the phase-shift A(j) that maximizes Eq. (3) also changes. For some masses, the optimal 
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phase-shift will occur at one of the discretely sampled time-shifts, and the results 
will be correct. For other masses, the true extrcmum will occur between discretely 
sampled time-shifts (as for 4096 Hz in the left panel), and the computed overlap will be 
erroneously too low. Henceforth in this paper, all overlaps are calculated at 32768 Hz. 
The LIGO/Virgo matched filter searches operate on data at 4096 Hz, however this 
issue is not a problem in these searches for reasons that can be seen in these two 
images. The loss of overlap by undersampling is no larger than 0.2%. The search 
utilizes a bank of templates which discrctizcs the parameter space. In constructing 
the template bank we have already incurred a potential loss of SNR of 3%, and this 
effect is therefore negligible. 

Before discussing overlaps between the submitted NINJA hybrid waveforms, we 
first present comparisons between time-domain post-Newtonian waveforms of the 
kind used to construct hybrid waveforms. (See the Appendix for a summary of the 
pN approximants used.) These waveforms terminate at Mw=0.136, 0.114, 0.069, 
and 0.135 for Tl, T2, T3 and T4 respectively At M = IOMq these correspond 
to termination frequencies of 435 Hz, 369 Hz, 222 Hz and 439 Hz, respectively. The 
results of our overlap calculations are shown in Fig. 6. Ref. [132] contains a far 
more extensive comparison between post-Newtonian approximants, although that 
analysis differs in several important aspects from what has been done here, notably the 
overlaps are maximized over the intrinsic mass parameters of one waveform. However 
the qualitative conclusion that Tl matches better with T2 and T4 than with T3 is 
consistent with our results. These results should be borne in mind when evaluating 
overlaps between hybrid waveforms using different pN approximants. In particular, 
overlaps at lower masses will be dominated by the influence of the two pN waverforms, 
and therefore hybrid waveforms using T3 will have relatively low matches against 
hybrid waveforms using Tl. Note also that T3 and T4 diverge from Tl as mass 
increases, this is precisely why we need to transition to numerical-relativity waveforms. 

Let us now discuss the main results of this section, the overlap between different 
submitted hybrid waveforms. For six of the black hole configurations listed in Table 1, 
hybrids have been constructed independently for different numerical waveforms, and 
overlaps were computed between each pair of waveforms in a group. Below we 
are showing overlap plots for all six cases: The q = {1,2} non-spinning waveform 
submissions are shown in Fig. 7, q = {3, 4} non-spinning waveforms in Fig. 8, and 
q = 1, Xi = 0.4,0.85 spinning waveforms in Fig. 9. We will discuss the q = {1,2} 
non-spinning cases in more detail. 

At the high-mass end, where the NR portion of the waveform dominates the 
overlap, the overlaps approach 1. Since all these waveforms model the same physical 
system, this is the expected behavior. High overlaps at high masses indicate good 
agreement between different NR codes, and the results here are consistent with the 
detailed study of the q = 1 case was performed in the Samurai project [21]. At the 
low-mass end, where the overlap is dominated by the pN portion of the waveform, the 
behavior is qualitatively as expected from Fig. 6. In particular the overlap between the 
BAM hybrid using Tl and the SpEC hybrid using T4 is lower than the matches between 
the other hybrids, all of which use T4. In the region between ^15 — 3OM0 there are 
drops in the matches between hybrids using T4. In this mass range the hybridization 
regions are passing through the most sensitive frequency band, and the mismatches 
are due in part to different choices of hybridization methods and parameters. The 
question of how many NR cycles are needed in order to produce a robust hybrid 
waveform is an area of active research [16-20]. 
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Figure 6. Overlaps between equal-mass, non-spinning, post-Newtonian time- 
domain waveforms. Note in particular the discrepancy between Tl and T4, as 
these were used in the majority of the hybrid waveforms. On the left overlaps 
arc calculated against the early aLIGO noise curve, on the right overlaps are 
calculated against the Zero-detuned, high-power noise curve. All waveforms 
include terms up to 3.5 pN-order, as do the post-Newtonian portions of the hybrid 
waveforms. 



If the approximants are the same, then the mismatches will depend only on the 
differences in (1) the hybridization methods, (2) the hybridization frequencies (and 
windows), and (3) the NR data. We have performed tests to verify that the overlap 
due to these effects is very small. We have made overlap calculations using a white- 
noise PSD, integrated between 10 Hz and 100 Hz, so the hybridisation region can pass 
fully into and out of band as the mass is varied between 10 Mq and 100 M©, and 
found that the maximum mismatch due to hybridization is 0.05% for the GATech 
and SpEC equal-mass nonspinning hybrids . This suggests that the contribution of 
the hybridization to the mismatch is very small, which is consistent with the results 
in [16]. Mismatches for masses M > 150 will be due only to differences between 
the numerical data, and we find these to be 0.1%, which is consistent with the results 
for equal-mass nonspinning binaries in [21] and for q = 2 nonspinning binaries in [124]. 

The overlap plots discussed thus far do not yet address the accuracy required for 
detection of gravitational waves. Broadly, in order to claim a detection a signal must 
match well against at least one template waveform used in a search. Several methods 
of assessing detection accuracy have been proposed [133, 134], however here we take 
the simple aproach that the waveforms arc sufficient for further detection studies if 
the overlap is above the the standard 0.97 threshold for a loss of no more than 10% of 
signals in a search. In some cases above the overlaps are not above this threshold, but 
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Figure 7. Differences between hybrid-waveforms for the same 
configuration. Plotted arc the overlaps of the SpEC (TaylorT4) hybrid against 
the four other submissions: All plots show results for both equal mass and 
mass ratio two hybrid waveforms with zero spin. On the top row overlaps were 
eomputed using the early aLIGO noise curve. On the bottom row the Zero- 
Detuned, High-Power noise curve was used. The overlaps are above 0.98 for 
identical pN-approximants. For different pN— approximants overlaps are above 
0.98 for the early noise curve and above 0.90 for the Zero-Detuned, High-Power 
noise curve. 



the appropriate quantity to evaluate to address this question is the overlap maximized 
over all of the physical parameters in a search. We now extend these overlap studies 
by maximizing over one physical parameter, the mass of one of the waveforms, as well 
as the time and phase. If the overlaps are now all over 0.97 (which they are), then 
they are acceptable for use in search-related studies. 

These extended overlaps also provide insight into the "parameter estimation 
question," as the error in parameter recovery between two hybrid waveforms gives 
a rough lower bound on the errors we may expect in recovering parameters of hybrid 
injections with search templates. In practice parameter recovery is likely to be worse 
than this, as it will involve maximizing over several parameters in addition to total 
mass, in addition to differences between the waveforms. Example plots using the 
equal-mass, non-spinning MayaKranc waveform as the signal and BAM hybridized with 
two different pN approximants as the template are shown in Fig. 10. Note that in 
the right panel of Fig. 10 the minimum overlap is below 0.97. This minimum match 
will be even worse for g > 1 binaries [16], but if maximization is done over the other 
physical parameters (mass ratio and spin) , then the overlap will increase to well above 
0.97 [20]. 

At the high-mass end the overlap is dominated by NR data, and as in Fig. 7 the 
overlaps are high without needing to move off the signal mass. At the low-mass end 
the same result would be expected in a pure pN/pN comparison although there is 
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Figure 8. Differences between hybrid-waveforms for the same 
configuration. Plotted arc the overlaps of the SpEC (TaylorT2) hybrid against 
the four other submissions: All plots show results for both mass ratio three and 
four nonspinning hybrid waveforms. On the top row overlaps were computed using 
the early aLIGO noise curve. On the bottom row the Zero-Detuned, High-Power 
noise curve was used. 



enough of the hybridization in-band to reduce the overlaps. However, changing the 
mass introduces a phase difference that accumulates over all the cycles in-band, and 
so higher overlaps cannot be achieved. The result is optimal mass values close to 
the correct mass value, but with a lower overlap. In the middle region these factors 
compete. At higher masses the overlap is reduced less by changing the mass and 
so the recovered value of the mass can stray further from the injected value. As 
the hybridization passes out of band this adjustment is no longer needed. The same 
general behavior can be seen in comparisons between non-spinning, unequal-mass 
{q = 2) waveforms, and for equal-mass spinning waveforms, as shown in Fig. 11 which 
shows overlaps between waveforms with identical parameters and hybridized with the 
same pN approximants maximized over the mass of one waveform. We can infer from 
these results that, although we have only made a crude estimate of the parameter bias, 
the accuracy of the waveforms (excluding the uncertainty in the pN approximants) is 
extremely high. There are many factors that may bias parameter estimation of real 
signals, including uncertainties in the detector calibration, noise in the detectors and 
errors in pN waveforms used as templates. Subsequent NINJA-2 data analysis studies 
will attempt to quantify the effects of these factors. 



6. Conclusion 



The efficiency of searches for GW signals from BH binaries and the accuracy with 
which the source parameters are estimated will depend crucially on progress in 
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Figure 9. Differences between hybrid-waveforms for the same 
configuration. Plotted arc the overlap comparisons of sets of equal mass 
spinning hybrids with Xi = 0-4, 0.85. On the top row overlaps were computed 
using the early aLIGO noise curve. On the bottom row the Zero-Detuned, 
High-Power noise curve was used. The overlaps are above 0.98 for identical pN- 
approximants. For different pN— approximants overlaps deteriorate to 0.94 for the 
early noise curve and above 0.85 for the Zero-Detuned, High-Power noise curve. 



incorporating information from approximation methods and numerical relativity into 
the waveform models used in data-analysis algorithms. Even in the nonspinning case, 
a recent extensive comparison of different pN approximants [132] concludes that for 
masses above 12 Mq numerical-relativity simulations of the last orbits and merger are 
required for the construction of optimal detection templates. This need is expected 
to be even greater when spinning binaries are included. And of course, numerical- 
relativity simulations will be yet more important for parameter estimation. 

Injections of inspiral-merger-ringdown (IMR) waveforms from analytical 
waveform models are already used to calibrate the analysis of detector data [135]. 
However, the synthesis of NR results into waveform models typically lags several years 
behind the most complete current sets of NR waveforms. The NINJA-2 project will 
consider direct injections of hybrid pN-NR waveforms, which will allow the use of 
waveforms which have not (yet) been used in constructing analytical models, and will 
avoid any additional modeling errors. Using the best available IMR waveforms will 
be particularly important for parameter estimation. 

The first goal of the second NINJA project has been to review the submitted 
hybrid waveforms. In this paper we have summarized the requirements that the 
submitted waveforms must meet (see Sec. 1). All submitted waveforms meet these 
requirements, within some minor caveats that are detailed in Sec. 1. We have also 
demonstrated the validity of the hybrid waveforms for use in the NINJA-2 project 
with more detailed consistency checks of the submitted data. 
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Figure 10. Overlaps and mass-bias when searching for one hybrid waveform with 
another hybrid waveform of the same configuration. Left panel: the g = 1 non- 
spinning MayaKranc waveform is taken as the signal, and the q = 1 non-spinning 
BAM waveform hybridized with TaylorT4 taken as the template. Right panel: 
the ? = 1, non-spinning MayaKranc waveform is taJjen as the signal and the g = 1 
non-spinning BAM waveform hybridized with TaylorT4 is taken as the template. 
In both panels, maximization is done over the mass of the template, as well as 
over time and phase. The horizontal axis gives the mass of the signal; the vertical 
axis gives the fractional difference between the injected mass and the mass of the 
template that maximizes the overlap. Overlaps are calculated against the early 
aLIGO noise curve. 



The validation of and comparisons between the submissions, which we have 
presented here, are a major feature of NINJA-2 that was not present in NINJA-1, 
and will be indispensable when analyzing the results from systematic injection studies 
over the course of the NINJA-2 project. A total of 63 waveforms from 8 different 
groups were contributed to NINJA-2, corresponding to 46 distinct NR waveforms, and 
29 different configurations of mass ratio and spins. For six configurations, multiple 
numcri(;al waveforms were submitted, and 16 numerical waveforms were hybridized 
with two different pN approximants (TaylorTl and TaylorT4). This has allowed 
waveform comparisons and tests of the accuracy standards, as discussed in detail 
in Sec. 5, and summarized below. For each submission, significant further information 
has been included in metadata files as specified in [60], e.g., hybridization frequencies, 
eccentricities, or literature references. This allowed us to automatically generate 
information such as Table 1 and has proved crucial in systematically analyzing the 
data set. Verifying these submissions, preparing a consistent data set, and evaluating 
their accuracy by comparing submissions with the same mass-ratio and spin, has been 
a formidable task, and in the present paper we only report on the {£,m) = (2,2) 
spherical-harmonic modes. Based on our various checks and comparisons as reported 
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Figure 11. Overlaps and mass-bias when searching for one hybrid waveform with 
another hybrid waveform of the same configuration. Left panel: the q = 2 non- 
spinning MayaKranc TaylorT4 waveform is taken as the signal, and the g = 2 non- 
spinning BAM waveform hybridized with TaylorT4 taJjen as the template. Right 
panel: the q = 1, Xi = X2 = 0.4 MayaKranc Taylor T4 waveform is taken as 
signal and searched for with the Llama waveform of the same parameters and pN 
approximant. In both panels, maximization is done over the mass of the template, 
as well as over time and phase. The horizontal axis gives the mass of the signal; 
the vertical axis gives the fractional difference between the injected mass and the 
mass of the template that maximizes the overlap. Overlaps are calculated against 
the early aLIGO noise curve. 



in this paper, we judge the submitted waveforms suitable for the NINJA-2 project. 

Parameter estimation places particular stringent demands on waveform templates. 
We believe the NINJA-2 waveforms to be of sufficient quality to estimate size and 
shape of maximum likelihood contours. However, the large truncation error of the 
post-Newtonian expansion in the cycles before and during hybridization will induce 
systematic errors; further investigations will be necessary before these waveforms can 
be used to check the accuracy of the estimated parameters themselfes (i.e. the location 
of the maximum likelihood contour, rather than its shape). 

As seen in Fig. 2, the submissions primarily cover only two lines in the {q, x) 
plane, leaving large regions of parameter space unexplored. We also do not consider 
any precessing signals, which adds several dimensions to the parameter space. These 
issues will be explored in future NINJA projects. Extending our analysis to non- 
dominant spherical-harmonic modes, and to a larger volume of parameter space, will 
mark a significant challenge for the future, and will require a significant extension of 
our methods of analysis, and of automatizing them, and will constitute a substantial 
research project by itself. 

An important shortcoming in our analysis of waveform accuracy is that it is 
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currently not possible to accurately estimate the pN truncation error. A reasonable 
approach seems to be to compare pN approximants, which are accurate to the 
highest known order, but differ in terms beyond this order (3.5pN without spins). 
These deviations still depend strongly on the choice of pN approximants used in any 
comparison, and on the mass ratio and spins. For example, in Figs. 6, 7 and 10 we 
compare the use of the TaylorTl and TaylorT4 approximants for mass ratio g = 1 
and zero spins, Fig. 7 also includes q ^ 1. Previous work implies that TaylorT4 
performs exceptionally well in these cases [29,35], apparently by coincidence, but 
beyond that it is not possible to make precise quantitative statements; this figure 
simply illustrates the general level of mismatch error that we may expect between 
various pN approximants across a range of binary masses. 

We caution that the comparisons presented in Sec. 5 were done only for subset of 
the black hole configurations, which do not include the most extreme cases, and not 
the cases with unequal Kerr parameters. Further work will be necessary to establish 
with similar confidence that all NINJA-2 hybrid waveforms are of similar quality. 

Our waveform comparisons in Sec. 5 are consistent with previous results: 
mismatches between waveforms are dominated by the choice of pN approximant, 
while NR and hybridization errors are far smaller. Hybridization choices and methods 
certainly affect overlaps, although the mismatch due to hybridization appears to be 
less than a fraction of a percent. This is not expected to have any noticeable effect 
on detection, but is likely to impact parameter estimation. The degree to which 
these effects will bias searches is one of the key question we hope NINJA-2 will be 
able to answer. In particular, theoretical studies such as the ones presented in this 
paper cannot replace a complete parameter estimation analysis, and it is not clear 
how theoretical studies based on Gaussian noise can predict the performance of GW 
searches in real noise. Another important goal of the NINJA-2 project is thus to build 
up experience in comparing simple theoretical studies with the full GW-search plus 
parameter estimation pipeline exercised in actual observations of GW events. 
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Appendix A. Post-Newtonian Waveforms 

The accuracy of hybrid waveforms depends very sensitively on the accuracy of the 
post-Newtonian waveforms with which they arc constriictcd. For this reason, the 
NINJA-2 collaboration invested significant effort into ensuring that all contributions 
used the most current pN information available. Mathematica notebooks containing 
the full expressions and derivations of the various approximants arc provided in the 
ancillary files available with this paper online. Here, we describe the techniques. 

The pN approximants used in this paper solve for the orbital motion to high 
accuracy by assuming that the motion is an adiabatic quasicircular inspiral, and by 
assuming that the loss of orbital binding energy E during inspiral is balanced by the 
flux of energy in gravitational radiation to infinity J" and the flux of energy going into 
the individual black holes caused by tidal heating M . The first assumption means that 
the motion can be described completely by the orbital phase function $(i). We then 
define the pN expansion parameter v := (Md$/dt)^/^, where M is the sum of the 
apparent-horizon masses of the black holes. The orbital binding energy, gravitational- 
wave flux, and tidal heating can be expressed as functions of this parameter. Thus, 
the energy-balance equation can be written ss E-\-J^-\-M = Q. Using the chain rule to 
rewrite E, we can rearrange this as an expression for dw/df, and include the expression 
for d$/dt to obtain a complete system of ordinary differential equations describing 
the motion of the binary: 

Av _ T{v) + M{v) d$ _ 

dt~ mv) d^~M' ^ 

The formulation of the pN approximants, then, comes down to writing down the 
expressions for E{v), J^{v), and M{v), and integrating the system for v{t) and ^{t). 
At leading order, the expressions for E and T are 

E{v) = -^'^ and ^(t;) = ft;i%^ (A.2) 

where higher-order terms include additional factors of v. The additional terms are 
currently known up to (3.5pN) for nonspinning systems, and to lower order for 
spinning systems. The tidal heating M is equivalent to a 2.5pN term in the flux 
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expression. The full expressions for E, J", and M are calculated in references [136- 
142], with a collection of errata in rcfc;r(^ncc^ [GO], and arc; given explicitly in the 
accompanying Mathematica notebook PNOrbitalPhase .nb, which can be found in 
the ancillary material with this paper online. For all the results in this paper, the full 
expressions were used. 

To integrate the system of ordinary differential equations, various methods have 
been developed, each of which should be equivalent at the level of accuracy of the 
pN approximation. These "approximants" have been given the names TaylorTl 
through TaylorT4 [29, 143, 144]. The TaylorTl approximant is constructed by directly 
evaluating the expressions in Eq. (A.l), which are then integrated numerically. For 
the TaylorT4 approximant, the right-hand side of the expression for Av/dt is first 
re-expanded in a Taylor series, and truncated at the appropriate 3.5pN order, which 
is then integrated numerically. The difference between these approximants is at the 
level of the 4pN uncertainty. 

Alternatively, we can find analytical formulas, rather than integrating 
numerically. The TaylorT2 approximant is derived by taking the (multiplicative) 
inverse of the first expression in Eq. (A.l) to construct a new system: 

dt _ E'{v) d$ _ d$ dt _ E'{v) 

dv ^ J^(v) + M{v) dv " dt dv~ M J^{v) + M{v) 

We re-expand the right-hand sides of these expressions in Taylor series, truncate at 
the appropriate 3.5pN order, and integrate with respect to v, obtaining expressions for 
t{v) and <I>(i'). This parametrically determines the phase of the binary as a function 
of time, with v being the independent parameter. Finally, the TaylorT3 approximant 
is constructed by inverting the series t{v) to obtain v{t). There is a 3pN logarithmic 
term in i(w), which must be treated as a constant in order to invert the series. Once 
this is done, the series for v{t) can be inserted into d$/di = v^/M, which can be 
integrated analytically to find ^{t). 

Using any of these approximants, the resulting orbital phase and frequency allow 
us to calculate the metric perturbation function h. The perturbation falls off at leading 
order as l/r and varies with angle, typically being decomposed in spin- weight s = —2 
spherical harmonics [60]. The dominant mode of this decomposition is generally the 
{£,m) = (2,2) mode. At leading order we have| 

/i2,2(t;,$) = — ,/f 87?t;2e-'2*, (A.4) 



r V 5 

where higher-order terms include additional factors of v. The additional terms arc 
currently known up to (3pN) for nonspinning systems [145] , and to lower order for 
spinning systems [60, 146]. The full expressions used for this paper, including other 
modes, arc given in the accompanying Mathematica notebook PNWavef orm.nb found 
in the ancillary materials with this paper online. 

f Note that most references discuss a change of variables given by 'I' := ^ — ln(v/vo), which is a 
relative 4pN modification to the phase (and can therefore be ignored at the level of accuracy currently 
known for the phase evolution). Here, vq is a freely specifiable parameter related to the origin of the 
time coordinate. This modification removes related terms in expressions for the waveform amplitude 
by shifting them to terms at 4.5pN order and higher, which are currently unknown. We emphasize 
that — at the level of current pN knowledge — this new variable ^ never needs to be calculated explicitly 
from the known orbital phase and frequency. Where the expressions for h refer to the standard 
orbital phase $ derived from the TaylorTn approximants may be used directly in place of * without 
subtracting the logarithm, and any term involving ln(ii/iio) should be removed from the expressions 
for amplitude. 
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The approximants just discussed describe the pN waveform in the time domain. 

For many purposes, it can be useful to have expressions in the frequency domain. 
The frequency-domain approximant TaylorF2 is obtained from TaylorT2 using the 
stationary-phase approximation (SPA) [147], together with the approximation that 
the gravitational- wave frequency is just twice the orbital frequency, so / = w^/ttM. 
Then the frequency- domain waveform is given by 



where i) is given in equation (A.l), and t{v) and ^{v) are the results of the TaylorT2 
approximant. 
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